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Improved lattice operators for non-relativistic fermions 
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In this work I apply a recently proposed improvement procedure, originally conceived to reduce 
finite lattice spacing effects in transfer matrices for dilute Fermi systems, to tuning operators for the 
calculation of observables. I construct, in particular, highly improved representations for the energy 
and the contact, as a first step in an improvement program for finite-temperature calculations. I 
illustrate the effects of improvement on those quantities with a ground-state lattice calculation at 
unitarity. 
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I. INTRODUCTION 



One of the main sources of uncertainty in Monte Carlo 
(MC) calculations of lattice field theories stems from fi- 
nite lattice spacing effects [1] . This is true in particular of 
non-relativistic systems at finite density, where the con- 
tinuum limit is approached by taking the limit of dilute- 
ness, such that the interparticle distance ~ kp^ (where 
kp is the Fermi momentum) is much larger than the lat- 
tice spacing I. This is achieved in practice by first con- 
sidering large volumes at fixed particle density, and then 
taking the zero-density limit. While this process is well 
defined, it is cumbersome to carry out, in part because 
calculations in large volumes require substantial amounts 
of CPU time. It is therefore useful to consider alternative 
approaches, based on effective field theory and renormal- 
ization group concepts [2 [3] , which treat lattice-spacing 
effects by modifying the ultraviolet (UV) dynamics of the 
theory. In such formulations, lattice-spacing effects are 
eliminated at finite volume, even for densities that are 
not low by conventional standards. 

In recent work, Endres et al. [31 [S] proposed a novel 
way to systematically reduce lattice-spacing effects in cal- 
culations of non-relativistic fermions. The method en- 
ables tuning of the lattice theory to high accuracy, such 
that the low end of the continuum energy spectrum is 
reproduced for the desired values of the two-body scat- 
tering parameters in the effective range expansion: 
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where S is the scattering phase shift, a is the scatter- 
ing length, and r^^ is the effective range. The connec- 
tion between the bare lattice theory (or rather its two- 
body spectrum) and these physical parameters is given 
by Liischer's formula [7 , which relates the phase shift to 
the energy E = p^ /m of the two-body problem in a box 
of side L: 
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where the sum is over all 3D integer vectors (the evalu- 
ation of S{ri) is discussed in Appendix A), and 9(x) is 
the Heaviside function. Throughout this work, we shall 
take units such that h — kg — m = 1, where m is the 
mass of the fermions. 

The two-body matching condition described above 
completely specifies the physics of dilute systems, such as 
those currently realized with ultracold atomic gases (see 
e.g. Ref. [^ for a review of the experimental situation). 
In this work we shall restrict ourselves to such systems, 
neglecting the additional complications that arise for in- 
stance at higher densities or for nuclear systems, where 
three- and higher-body forces play an important role. 

We shall briefly review the work of Ref. |3] below, 
but it is useful to mention the main steps underlying 
the method at this point. First, one writes down the 
transfer matrix T, representing the two-body interaction 
via a generalized Hubbard-Stratonovich (HS) transfor- 
mation [SI [5]. The latter contains N^rj arbitrary coef- 
ficients C„. The matrix elements of T are then com- 
puted in the subspace of two particles at vanishing total 
momentum. The resulting matrix is then diagonalized 
in the s-wave subspace (or rather its lattice equivalent), 
and the HS coefficients C„ are tuned using an iterative 
algorithm such that the eigenvalues of the proposed T 
match exp(— r_E), where r is the temporal discretization 
parameter and E are the energies dictated by Liischer's 
formula, for the desired box size L and choice of scatter- 
ing parameters. 

In the case of the unitary limit, where by definition 



p cot S{p) = 0, 
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the above procedure yields considerable improvement in 
the approach to the continuum. Indeed, with a single co- 
efficient Cg (the simplest case, discussed in Sec. In]) one 
may tune the scattering length a, but the effective range 
r^g remains finite due to lattice-spacing artifacts. Natu- 
rally, by introducing more parameters one may tune the 



effective-range expansion in Eq. (11) to higher accuracy, 
thus systematically eliminating the need for extrapola- 
tions to the dilute limit and leaving only finite-volume 
effects unaccounted for. 

It is the main objective of this work to extend the ap- 
proach of Ref . 4 to designing not only improved transfer 
matrices but also improved operators for the calculation 
of observables. This is an important step toward reducing 
lattice-spacing effects in finite-temperature calculations. 
While these effects were studied at unitarity in Ref. [10] 
for the Hubbard model using dynamical mean-field the- 
ory, careful systematic studies in full-fledged MC calcu- 
lations have only recently started to appear [21 [3 [TTl [T^] 
and are still restricted to the ground state. Interestingly, 
recent studies based on functional renormalization group 
techniques [13, have also started to tackle the problem of 
extrapolating to infinite volume and particle number in 
the unitary limit. 

In Sec. ITT] we analyze a basic example, as a primer to 
reviewing the more sophisticated formalism of Ref. [1], 
which we explain briefly in Sec. 



Ill 



In Sec. IV we present 
the main developments of this work, with the corre- 
sponding illustrative results and conclusions appearing 
in Sec. IV] Finally, some of the less trivial numerical is- 
sues are explained in the Appendix. 



II. A SIMPLE EXAMPLE 

Before proceeding to a more detailed discussion, let us 
analyze a simple case, to fix notation as well as ideas. 
Consider the following lattice Hamiltonian: 
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where s denotes the spin projection, g is the bare lattice 
coupling constant, and n^ ^ = a], ^a^ ^ denotes the num- 
ber density operator at lattice position i for spin projec- 
tion s. The value of g is tuned to the desired two-body 
scattering properties by solving the two-body problem on 
the lattice and finding the scattering amplitude, which in 
this case can be done analytically. 

The transfer matrix is then expressed approximately 
in powers of the imaginary time step r using a Suzuki- 
Trotter decomposition, for instance as follows: 
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The kinetic energy factor in Eq. ^ is easy to treat, as it 
is an exponential of a one-body operator which is diago- 
nal in momentum space. Notice that we have chosen to 



define the dispersion relation E — k'^ /2m in momentum 
space, as opposed to defining it via a discrete representa- 
tion of the Laplacian operator in coordinate space (which 
is common in Hubbard- model type formulations). 

On the other hand, the potential energy factor repre- 
sents a challenge, as it is the exponential of a non-trivial 
two-body operator (as the original Hamiltonian). It is 
at this point that the HS transformation plays a central 
role, allowing us to exchange the complexity of the two- 
body operator for a path integral involving only one-body 
operators. Specifically, we write 

exp(-T¥^)J|exp(r5n^ -fi^.) (9) 
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where A = 2{e'^^ — 1), Va — Y[idai/{2TT), and (Ji is an 
external auxiliary field. In this way, one decouples the 
transfer matrix for each spin, such that we may write 

(10) 
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which in the one-particle subspace reduces to 
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Notice that in Eq. ^ we have used a version of the 
HS transform due to Lee [19] , in which the auxiliary field 
is continuous and compact (the integral is restricted to 
the interval [— tt, tt]), as opposed to more conventional 
versions that are discrete, or continuous but unbounded. 
Apart from the path integral, at this point applying the 
transfer matrix becomes a problem of applying a product 
of one-body operators, which one can easily deal with. 

All of the above steps are standard in the literature of 
many-body MC calculations, except perhaps for the use 
of a "perfect" dispersion relation defined in momentum 
space, which has become more common only in recent 
years [TT] [T5] . This feature can be regarded as a basic 
kind of improvement, as it reduces lattice-spacing effects 
relative to Hubbard-model approaches, where dE/dk — > 
at high momenta. 

This example captures the main features of most mod- 
ern many-fermion MC calculations. Nevertheless, the 
simplicity of this approach is in some ways excessive, 
largely because we only have one coupling constant at 
our disposal. Indeed, should we desire to fix more than 
one coefficient in the effective-range expansion, we would 
need a richer bare interaction, and a correspondingly 
more sophisticated HS transformation. In this simple 
case, to reach for instance the unitary limit, we can tune 



the scattering length, but we also need to consider very 
dilute systems to avoid the effects of finite range induced 
by the UV lattice cutoff it /I, as in that case one has 
r^ff ~ 0.40Z. Moreover, taking for example 80 particles in 
a 10"^ lattice volume, which shall be our illustrative ex- 
ample below, one finds kpr^g ~ 0.54, which is not nearly 
as small as one would like. 

As mentioned before, one can resort to alternative 
approaches that modify the Hamiltonian by including 
higher-order terms in a low-momentum expansion, tun- 
ing the corresponding couplings to eliminate UV effects. 
The latter strategy is essentially what Ref. [3] advocates, 
following the spirit of the Lattice QCD program initiated 
by Symanzik many years ago [3] to design improved effec- 
tive actions that approach the continuum limit (see e.g. 
Ref. 14 ). Since it will be useful for the rest of this work, 
we shall briefly review the method of Ref. [4 in the next 
section. 



III. IMPROVED TRANSFER MATRICES 

The work of Ref. [1] tackles the problem of reducing 
UV lattice effects by defining an improved transfer ma- 
trix, using a generalized HS transformation and a low- 
momentum expansion. This is accomplished by allowing 
the constant A to have a non-trivial momentum depen- 
dence. For simplicity, it is useful to take A{p) to be diag- 
onal in momentum space. Following Ref. [4J, we expand 
A{p) using a set of operators as 
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where one may choose to define, for convenience, 

0„(p)=fl-e-P'^V\ (14) 



as done in Ref. [3] (up to an n-dependent constant in 
front), or, as we shall use in the rest of this work. 



0„(p) = [2sin(p/2)] = 



(15) 



Notice that this last expression contains only even pow- 
ers of p = vp^ , and is therefore an analytic function of 
the momentum p. Both of the above choices for C'„(p) 
behave as ~ p^" at low momenta. For the purposes of 
Monte Carlo calculations, the operator ^(p) can be com- 
puted once and for all at the beginning of the calculation, 
and is applied at each time slice using Fourier transforms. 
In order to determine the coefficients C„, we find the 
explicit form of T for two particles, which upon integrat- 
ing the auxiliary field is given in momentum space by 
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where V is the lattice volume and T{k) = (fe2 + fcf)/2m. 
One cannot fail to notice that the above transfer matrix is 
not Galilean invariant. We elaborate on this issue in the 
appendices. In this work we shall take A = 7r(l — 10~^), 
following the steps of Ref. [4j. The operators in Eqs. (14) 



and (151 are defined to be constant above p = A , and 



the kinetic energy factor exp(— Tr/2) is defined to vanish 
above that boundary, such that there is no propagation 
for p^ > A^. 

Evaluating 7^ in the center-of-mass frame, we have 
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where p^. and q^ are incoming and outgoing relative mo- 
menta. By diagonalizing this expression, we may identify 
the eigenvalues of 7^ with e~'^^ , where E are the two- 
particle eigenvalues of the Hamiltonian we are implicitly 
defining. 

One may then tune the C„ such that this energy spec- 
trum matches the one required by Liischer's formula 
for the lowest Ni^ eigenvalues, given the desired val- 
ues of the scattering parameters in Eq. (fl]), and using 
E — p^ /ni — rf{2TT)'^/{mL'^). Specifically, if we are inter- 
ested in describing the unitary limit, the eigenvalues are 
determined by the zeros of the function S{r]) in Eq. (|3|. 
The actual fitting of the C„ can be performed iteratively, 
as described in detail in Ref. [3]. Notice, in particular, 
that our expression for 7^ ceases to be Hermitian when 
we promote A to be an operator. Therefore, one needs to 
use 72^7^ rather than 7^ to diagonalize and fit the eigen- 
values, as well as for the actual Monte Carlo calculations. 

The results of our fits for C„ are shown in Table [l] The 
quality of the improvement can be assessed by plotting 
pcot(5(p) as a function of rj'^, which is shown in Fig. [u 
for N^ = 20, r = 0.05, and Nq = 1 - 5. As can be 
appreciated in the figure, the same effect is achieved as 
in Ref. |3]: as the order of the expansion is increased, 
the transfer matrix is accurately tuned to unitarity up to 
progressively higher momenta. 



IV. IMPROVED OBSERVABLES 

The above procedure represents a significant step for- 
ward in mitigating lattice-spacing effects in MC calcula- 



TABLE I. Results of fitting the coefficients C„ (see Sec. Ill I 
to the low-energy spectrum of the two-body problem at res- 
onance, in a box of side N^ = 16, for an imaginary time step 
T — 0.05, in lattice units. 
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FIG. 1. (Color online) Plot oi pcot5{p), in lattice units, as 
a function of ri^ = -B2/P0 (where po = 27r/L), for A^'^ = 20, 
r — 0.05, and levels of improvement Nq = 1 — 5, for a Galilean 
non-invariant definition of the transfer matrix (for Nq > 2) . 



tions, especially considering that it requires only a small 
coding investment for its implementation in extant MC 
codes, and it results in minimal computational overhead. 

A somewhat unsettling issue remains, however, partic- 
ularly in connection with improving finite temperature 
lattice calculations, such as those of Refs. [T5H17) . In- 
deed, in those calculations, as well as in similar ground- 
state approaches, the transfer matrix is not the only ob- 
ject carrying lattice-spacing effects: the operators used 
to compute expectation values also suffer from the same 
problems. The situation may appear problematic at first 
sight, as the improvement strategy defined above does 
not directly define improved operators that we could use 
to compute observables at finite temperature. 

On the other hand, in conventional formulations, 
knowledge of the explicit r dependence of the transfer 
matrix (as in the simple example described in Sec. In]) al- 
lows us to take a derivative that brings the Hamiltonian 
down from the exponent (c.f. Eqs. (10 1 and (11)), which 



results in a practical expression for the calculation of the 
energy [? ]. 

To fix ideas, let us consider the grand canonical ensem- 
ble (see also Ref. [El [18]), where the partition function 
satisfies, by definition. 



-pn 
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Averages of observables can be extracted from knowledge 
of Z through derivatives, and in particular the energy is 
obtained by means of 

dlogZ 

w 

When using the formalism derived in Sec. |llj we have 
(assuming the system is unpolarized) 



= E- jjN. 
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Z = Tr [r^-] = jVa dct [1 + U[a] 
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and where now a is to be regarded as a space-time varying 
field, and cr^ is a restricted to the t-th imaginary-time 
slice. The determinant is to be taken in the subspace of 
one-particle states, and 



U[a]^ll%.[<j, 



(23) 



As the derivative in Eq. ( |19[ ) becomes now a derivative 
with respect to t, it is clear that all we require is knowing 
how to differentiate 7^ with respect to r, because 
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We shall take the point of view that differentiation 
of partition functions with respect to parameters in the 
Hamiltonian is the proper way to obtain expectation val- 
ues of operators. We shall then generalize the improve- 
ment procedure of Ref. W to design highly improved op- 
erators (in particular for the energy and the contact) and 
accomplish this by taking formal derivatives of the trans- 
fer matrix and writing them in an operator expansion, fit- 
ting the expansion coefficients to reproduce the low end 
of the exact two-particle spectra. By " highly improved" 
we mean improvements that go well beyond tuning the 
first two coefficients in the effective-range expansion. 

While in the above example we have focused on the 
calculation of the energy, which we resume in the next 
section, similar derivations apply for the calculation of 
the contact, as we shall see in Sec. |IVB| 



A. Energy 

In the simple case presented in Sec. |TT1 the calculation 
of the first T derivative results in the following expression: 
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where K is the following anticommutator 



and 
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The derivative of \/ A with respect to r appearing in 
Eq. ( 28 ) is known analytically in this case. When using 
improved transfer matrices, however, that derivative is 
somewhat harder to compute, as we only know the co- 
efficients C„ in Eq. (13 1 as a result of a numerical fit. 



which complicates the calculation of dCn/dr. 

Here we extend the procedure of Endres et al. to fitting 
the coefficients in the expansion of the derivative of ^(p). 
There are at least two advantages in following this route 
over computing the derivatives via finite differences. In 
the first place, the fitting procedure can provide a much 
more accurate determination of the expansion coefficients 
of the derivative, with less effort. (That being said, fi- 
nite differences do provide a useful starting point for the 
fitting routine.) Secondly, this route allows one to use 
different operators for different observables, regardless of 
what is used for ^(p) itself. This feature can potentially 
be very convenient: different observables are in general 
sensitive to different regions of momentum space, such 
that there is no a priori reason why a given set of opera- 
tors On should be universally useful. 

The main idea behind the method remains the same as 
in Ref. [Ij. Liischer's formula provides us with the exact 
two-particle spectrum, such that the eigenvalues of the 
exact two-particle transfer matrix and its r derivative (s) 
are known; 
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where £'2 ^^^ the exact two-particle energies in a contin- 
uous box [? ]. 

In order to match this spectrum, we take the derivative 



of Eq. (17 1 with respect to r, evaluated between eigen- 
states l-E) of the proposed transfer matrix 7^. We thus 
obtain, using the Feynman-Hellmann theorem, 
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The rest of the recipe consists in taking the right-hand 
side of Eq. ( 30 ) and fitting the coefficients Z3„ such that 



the first A^„jax eigenvalues of the exact expression Eq. ( 29 ) 



(which correspond to the lowest eigenvalues prescribed by 
Liischer's formula) are reproduced. We may assume at 
this point that the coefficients C„ are known, such that 
the eigenvectors \E) are fixed when we set out to find 
the Z?„. In that case, the D^ are determined by a linear 
system of equations of order iVmax x -^max^ 
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TABLE II. Results of fitting the coefficients D,i (see Sec. IV I 
to the low-energy spectrum of the two-body problem at res- 
onance, in a box of side N^ = 16, for an imaginary time step 
T — 0.05, in lattice units. 
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FIG. 2. (Color online) Logarithmic plot of the difference 
between the exact spectrum of the two-body problem Eq. ( 29 1 
at unitarity, and the approximate spectrum obtained with 
various levels of improvement A''^ = 1 — 5, relative to the exact 
spectrum (see text for details), as a function of 77^ = -E2/P01 
where po = 2n/L. The original data are discrete; the lines 
are intended as a guide to the eyes. These results correspond 
to a lattice of side Nx ~ 20 and temporal spacing r — 0.05. 



where 



and 



MEn^^{E\OjE) 
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YE^Eexp{~TE)~{E\e-^K2e-^\E) (35) 

Once the coefficients Z?„ have been determined, one 
can use Eq. ( |28[ ) in a lattice calculation simply by replac- 
ing A — > A(pJ7and taking dCj^/dr = D„. The results of 
the fits for C„, D„ are shown in Tables |I] and |TT] 

As a first illustration of the level of improvement that 
can be achieved for the energy, Fig.[2]shows the difference 
between the approximate spectrum with various levels of 
improvement i^approx and the exact spectrum £^exact, as 
a function of 77^, through the quantity Log(|A£'|), where 
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t)/Ee- 



As expected, with each 



new parameter a new eigenvalue is reproduced, with the 
concomitant reduction in the error. 



Surprisingly, the approximations perform well consid- 
erably beyond the lowest eigenvalues that are fit. Indeed, 
with each new level of improvement we see, apart from 
a dramatic reduction in the error for the target eigenval- 
ues, an extra reduction in the error that is evident even 
beyond 77^ ~ 30. 



B. Contact 

One of the many ways to define the contact C (see 
Refs. [20I [21]) is through the derivative of the energy 
with respect to the inverse scattering length: 
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Since E may be obtained directly from the logarithm of 
the partition function, we are again in a situation where 
we require a derivative of the transfer matrix with respect 
to a parameter, in this case a~^. The generalization of 
the above definition to finite temperature in the grand 
canonical ensemble is 



9^ 
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where Vl is the grand thermodynamic potential, T is the 
temperature and /x is the chemical potential. 

In many-body lattice calculations, using these defini- 
tions involves the following expression: 
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The first step towards using these expressions in com- 
bination with the improvement procedure is to take a 
formal derivative of T with respect to a^^ in the two- 
particle space, which we can treat exactly: 



n'T-cxact 

aa-1 



T^exp(-rii;2). (40) 



In order to compute the change in the exact two-particle 
energy E2 due to a small change in the inverse scattering 
length, we use the fact that the energies are implicitly 
defined as solutions of Eq. ([2| , which implies 



dE^ 






dS 
drf- 



(41) 



where 77^ = i?2i^/(27r)^, and the derivative on the right- 
hand side is to be evaluated at the corresponding solution 
of Eq. ([2]). Table [lV| in the Appendix shows the first 30 
roots of S{iiii) and the corresponding values of dS/drj-. 



TABLE III. Results of fitting the coefficients F„ (see Sec.jIVl 
to the low-energy spectrum of the two-body problem at res- 
onance, in a box of side A'^^, = 16, for an imaginary time step 
r — 0.05, in lattice units. 
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-0.00023 



0.00180 



In the derivation of Eq. (41 1, we have assumed that all 



the effective-range parameters other than the scattering 
length are kept constant. 

Having the exact target spectrum, we proceed by find- 
ing the corresponding expression in terms of the HS 
function v4(p), which we obtain using Eq. (17 1 and the 
Feynman-Hellmann theorem: 
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where, as before, we expand in terms of our chosen set of 
operators. 



dA{p,) 
5a-i 



N„ 



E FnOniPr) 



(43) 



and we determine the coefficients F,^ by fitting the diag- 
onal matrix elements in the right-hand side of Eq. ( 42 ) 
to the exact spectrum of Eqs. (40) and (41). As with the 



energy, the fitting procedure can be reduced to solving a 
set of linear equations of order iVmax x A^max- Illustrative 
results of such a fit are shown in Table IIIII 

As with the energy, a first glimpse at the level of im- 
provement that can be obtained at this point. This is 
shown in Fig. |3] where we display the difference AC be- 
tween the spectrum with various levels of improvement 
and the exact spectrum, divided by the latter. As in 
Fig. [21 each new parameter allows one to fit a new eigen- 
value to high accuracy, matching the desired physics be- 
yond the lowest momentum. 

Unlike in Fig. [2J the improvement for eigenvalues be- 
yond those explicitly fit is limited, breaking down after 
the eighth or ninth eigenvalue. From that point on to- 
ward higher energies, little improvement, if any, is ob- 
served as the order of the expansion is increased. This 
behavior is only unexpected in the light of Fig. [2] where 
the situation is (surprisingly) much more favorable. Pos- 
sibly, part of the reason for this behavior is that the tar- 
get spectra for the energy and the transfer matrix are 
monotonic (either increasing or decreasing) at low ener- 
gies, whereas the one for the contact is closer to a random 



sequence (see Table IV in the Appendix). 

On the other hand, as mentioned before, there is no 
reason to believe that a given choice of operators will be 
equally useful for all observables. It is therefore natural 
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FIG. 3. (Color online) Logarithmic plot of t he difference 
between the exact two-body spectrum Eq. ( 40 1 at unitarity, 



and the spectrum obtained with various levels of improvement 
Nq = 1 — 5, relative to the exact spectrum (see text for 
details), as a function of rj^ = E2/P0, where po = 27r/I/. The 
original data are discrete; the lines are intended as a guide to 
the eyes. These results correspond to a lattice of side A^^ = 16 
and temporal spacing r = 0.05. 



to suspect that a different set of operators could yield a 
better expansion for the contact. This possibility remains 
to be studied. 



C. Extrapolation to the ground state 

Taking a Slater-determinant state \ipo) as a starting 
point, it is easy to see that the probability sum in ground- 
state calculations, for a temporal extent /3 can be written 
as 



Zo(/3) = ^Afee- 



-/3S, 



where E^. are the exact energy eigenvalues and 

A,^|(Vo|i?fe)|'. 



(44) 



(45) 



Taking a derivative of logZo(/3) with respect to r it is 
easy to see that in the large-/? limit one can write 



EW) 



dp 



^En + b^e-^\ 



(46) 



where Eq is the ground-state energy, S ^ Ei — Eq, and 



6^ = 



A, 



-iE,-E,). 



(47) 



Similarly, taking a derivative with respect to a~^ one 
can see that the first few leading contributions to the 
asymptotic behavior at large /3 are given by 

CiP) ^ |f^^l|# - C„-H6,,r^+6,,e-^^- (48) 



where Co is the ground-state contact, 
Atttti d log Aq 

Airm A^ ( dE^ 



bci = 



bc2 = -- 



dEn 



h^ An V9a-i g^-i 



(49) 
(50) 



We shall use these expressions to motivate the extrapo- 
lations to the /3 — >■ 00 limit below. 



V. ILLUSTRATIVE RESULTS AND 
CONCLUSIONS 

To illustrate the effect of improved operators in real- 
istic lattice calculations, this section presents the results 
of ground-state MC calculations of the energy (Fig. Ill) 
and the contact (Fig. ^, for an unpolarized system at 
unitarity. 

Results are shown for 80 particles (40 per spin) in a 
volume of 10'^ lattice points, for various levels of im- 
provement Nf^ = 1 — 4. In each case, calculations were 
performed for time directions of extent Pep = 2.0 — 8.0 
(corresponding to 7V^ roughly between 40 and 200), sub- 
sequently extrapolating to the /3 — >■ 00 limit. We have 
taken r = 0.05 in lattice units, and a Slater determinant 
of plane waves as the starting guess for the ground-state 
wavefunction. For each value of /3, we obtained approxi- 
mately 400 samples of the auxiliary field a. The statistics 
is enhanced by a factor of 20 — 100 by the fact that the 
operators we have defined utilize every other time slice. 
Some obvious fluctuations remain, as evident from some 
degree of jaggedness in the data. This can be resolved by 
increasing the statistics, but it does not affect our main 
conclusions. 
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FIG. 4. (Color online) Energy, in units of EpQ — ^Nep, as 
a function of the extent of the time direction fiep = tNt^p, 
for 80 particles and levels of improvement Nq = 1 — 4. 



Once the improvements are turned on, i.e. for No = 
2 — 4, we see that the change from kpr^g ~ 0.54 to 
kpT^a = results in a considerable reduction in the 
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FIG. 5. (Color online) Contact, in units of Nkp, where 
kp = {Sn'^N/Vy''^ is the Fermi momentum, as a function of 
the extent of the time direction /3e^ — rN-rep, for 80 particles 
and levels of improvement Nq = 1 — 4. Also shown are fits 
to the Nq = 1 and Nq — 4 datasets using the first two terms 
in the asymptotic form of the previous section (see text for 
details). 



energy. Our extrapolations to the ground state yield 
^ = E/EpQ = 0.467(2) for the unimproved case, and 
^ = 0.402(1), 0.399(2), 0.399(2) for No = 2,3,4, respec- 
tively. This change, of roughly 15%, is in line with 
our expectations based on the large-scale calculations of 
Ref. [TT] . While these results are still roughly 5% greater 
than those of Ref. [IJ, which found ^ « 0.38, it is re- 
markable that this can be achieved with a 10^ lattice. 
The remaining effects must be a combination of finite 
volume and finite imaginary-time step r effects. 

A similar situation is observed for the contact. Using 
the extrapolation formula of the previous section at lead- 
ing plus next-to leading order, we obtain C/{Nkp) = 
3.17(3) in the unimproved case, and C/{Nkp) = 
3.31(4), 3.34(4) and 3.30(4) for No = 2, 3, 4, respectively. 
(We have discarded three data points at the lowest val- 
ues of f3ep for the purpose of capturing the asymptotic 
behavior at large /3; the fits are stable against further 
removal of points at low /3ep). While the latter extrapo- 
lations clearly overlap when taking the uncertainty into 
account, they do not overlap with the unimproved case, 
which is clearly consistent with what we see in Fig. [5] at 
large /3. The remaining systematic errors, likely due to 
volume effects for the most part, appear to be only as 
large as 3%, if we consider the most recent estimates of 
the contact in the ground state (« 3.39, see Ref. |24|). 
As in the case of the energy, it is remarkable that such a 
small volume as 10^ already yields a result that is quite 
close to the best current estimate. 

In conclusion, we have presented a methodology based 
on the work of Refs. [11 [5], whereby one can generate not 
only improved transfer matrices but also improved opera- 
tors, which account for finite-range effects in a systematic 
fashion. The improvement program can be carried out in 



a Galilean invariant or in a Galilean non-invariant way. 
We have chosen the latter because it provides a better 
approach to the unitary limit. As shown in Appendix [Bj 
the difference in the eigenvalues of the invariant and non- 
invariant transfer matrices is very small, such that the 
effect can be safely considered to be negligible. This 
is likely due to the fact that the lattice already breaks 
Galilean invariance. With the chosen improvements in 
place, we see a large change for the energy but a rather 
small change for the contact. Indeed, once finite-range 
effects are under control, the dominant contribution to 
systematic uncertainties is expected to be given by finite- 
volume effects. If so, the latter appear to be relatively 
small both for the energy (roughly 5%) as well as for the 
contact (roughly 3%). 

While we have focused on the unitary limit, the tun- 
ing procedure for the transfer matrix and the operators 
can easily be applied to systems away from unitarity as 
well as away from the zero effective-range limit. In com- 
bination with hybrid Monte Carlo (HMC) techniques, 
which have recently made it possible to access larger 
volumes than ever before (as large as 20^), we expect 
the method presented here to provide a powerful strat- 
egy to tackle the non-relativistic many-fermion problem. 
This applies in particular to finite-temperature calcula- 
tions where lowering the ultraviolet cutoff (while keeping 
effective-range effects under control) reduces the size of 
the basis, speeding up the computations considerably. 

In this regard, it should also be pointed out that im- 
proving the transfer matrix affects the performance of the 
HMC algorithm only marginally. Indeed, the scaling with 
system size (particle number, spacetime volume) remains 
unchanged when increasing No- Only the prefactor in 
the scaling law increases somewhat (by about 20 — 25%) 
due to the need for extra Fourier transforms when apply- 
ing the operator A{p). Aside from this, the force in the 
molecular dynamics (MD) part of the HMC algorithm 
becomes somewhat larger when A(p) is improved, such 
that somewhat smaller MD time steps (not to be con- 
fused with r) are needed. In short, the performance of 
the HMC algorithm is largely unaffected by the use of 
improved transfer matrices. 

Finally, we would like to stress that the operator pro- 
posed here for the calculation of the contact (based on 
the derivative of the two-body energy, in turn computed 
using Eq. (41)) constitutes a novel way to calculate C in 
MC calculations. This method should be contrasted with 
others based on extracting the derivative of the equation 
of state via finite differences, using the momentum distri- 
bution, or using the pair distribution function. Each of 
these will, in general, behave differently in terms of their 
systematic effects. 
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Appendix A: Evaluation of 5(r;). 

There are many ways to evaluate the function S{ri) ef- 
ficiently. The one that I have found easiest to implement 
and understand is essentially the one communicated to 
me by Shina Tan, which I reproduce here, with some 
minor modifications that I have found useful. 

The first step is to enhance the convergence of the sum 
by introducing an exponential factor and separating the 
singularity in the tail (see Ref. [22]): 
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|n|<A 



n^ _ j^ 



E 

|n|<A 



1 



-xin'-rf) 



n? — if' 



E 

|n|<A 



-x{n 



^2 _ ^2 



(Al) 

where n = |n|, the sums are over all triplets of integers 
n, and a; is a small positive parameter. The second sum 
has no convergence problems, in fact it converges very 
quickly, so we shall put it aside by defining 



TABLE IV. 
tliose roots. 



First 30 roots of 5(r;), and dS/drf evaluated at 



9 

10 
11 

12 
13 
14 
15 
16 
17 
18 
19 
20 
21 
22 
23 
24 
25 
26 
27 
28 
29 
30 



~nk 



-0.0959007 
0.4728943 
1.4415913 
2.6270076 
3.5366199 
4.2517060 
5.5377008 
7.1962632 
8.2879537 
9.5345314 
10.5505341 
11.7014957 
12.3102392 
13.3831152 
15.3537375 
16.1218253 
17.5325415 
18.6053932 
19.5186394 
20.4033187 
21.6944179 
23.0194727 
24.3306210 
25.3016129 
26.6803600 
27.8780019 
29.6156511 
31.3536974 
32.1958982 
33.4483351 



dS/d4, 



123.82387 
39.75514 
82.36519 
106.24712 
84.23133 
161.88763 
212.49220 
62.95336 
231.79580 
247.82611 
233.82976 
185.61411 
183.65019 
316.68684 
82.86757 
506.59914 
371.40245 
308.00372 
255.97969 
329.98905 
394.81924 
94.98929 
342.25749 
526.27127 
514.90705 
150.20773 
548.38017 
114.02114 
443.21169 
452.78989 
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The first sum, on the other hand, converges just as 
slowly as the original one (once we subtract the 47rA 
term). Let us then analyze this first sum. For < a; <C 1, 
we may approximate it very accurately in terms of an in- 
tegral: 



where the 47rA cancels with the — 47rA term in Eq. (pi, 
and the integral in the second term is equal to 1/2^Jtt/x 
in the limit A — )• oo. Second, we have 



l = An dn- 
Jo 



11? — rp 



To treat this term we write 
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where we have used the spherical symmetry of the inte- 
grand. Writing 



T 



V? — rf 



9 9 ' 



(A4) 



we separate the integrals into two terms. First, 

/•A pA 

47r dn(l- e-^^"'-"')) = 47rA - 4^^"'/ dn e" 



(A5) 



1 _ e-^("'-''') = a;(n2 - rj^) / dy e-^^^"'""') (A7) 

and exchange the order of integration of y and n, to ob- 
tain 



pi pA 

I = Attxtj^ / dy e^^''' / dn e~"=^"' 
Jo Jo 



S{7j) =S{7],X)^ 



27r3/2 



e^"' - 2xrj^ / di e^" * , (A9) 



which is to be evaluated for < x ^ 1. 



(A8) 



all of which can be easily evaluated in the limit A — > oo, 
where the full result can be written as 
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With the above expressions at hand, it is straightfor- 
ward to find the first derivative of S with respect to i]^ : 



dS _ dS 
dif' diif' 



where 



(AlO) 



dS ^, , , V- e^^t"^" ) 

' |n|<A ^ ' ' 



(All) 



Table IV shows the first 30 roots of 5(?y), along with 
the corresponding values of the derivative dS /drj^ . 



Appendix B: Galilean invariance vs. non-invariance. 



One cannot fail to notice that not only is the transfer 



matrix of Eq. (16) not symmetric, an issue we managed 



to deal with above, it is also not fully Galilean invari- 
ant. Indeed, while translation and rotation invariance 
are preserved (in their discrete forms set by the lattice), 
Galilean-boost symmetry is broken. This is because we 
have assumed that A was promoted to an operator that 
acts on everything that appears to its right. 

On the other hand, if we assume that A acts on the 
auxiliary field only, as in Ref. [J, one arrives directly 
at a local, Galilean-invariant lattice theory. Indeed, in 
that case Eq. ( 16 ) changes form in a very simple manner, 
namely the function A{pA is replaced by A{p^—q^), and 
similarly for the other spin. The function A has then a 
clear physical significance as the momentum transfer in 
a two-body collision. 
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FIG. 6. (Color online) Plot of pcot5(p), in lattice units, as 
a function of rf = -E'2/po (where Po = 27r/L), for Nx = 20, 
r = 0.05, and levels of improvement Nq = 1 — 5, for a Galilean 
invariant definition of the transfer matrix. 



Apart from the obvious desirability of Galilean invari- 
ance, these observations seem to indicate that one should 



choose this form over the non-invariant version. How- 
ever, when comparing the approach to the unitary point 
based on improvements, as shown in Figs. [T] and [6j it is 
clear that the non-invariant version is much superior to 
its invariant counterpart. Furthermore, comparing with 
the results of Ref. |4] , it seems clear that one would need 
much larger lattices in the Galilean invariant version to 
make Fig. |6] look like Fig. [T] which we have checked with 
our codes. (Note that Fig. 2 in Ref. [1] presents the same 
kind of plot but corresponding to a much larger lattice 
size, namely N^ = 32.) 

The reason for this difference is likely due to the fact 
that (discrete) Galilean invariance is only respected in an 
infinite lattice volume. Indeed, in any finite lattice the 
center-of-mass and relative motions are actually not sep- 
arable. Boosting to frames of different total momentum 
will, in general, result in Hilbert spaces of vastly different 
dimensions for the subspace of relative motion. In par- 
ticular, the zero total momentum Hilbert space, used to 
tune our two-body interaction, has by far the highest di- 
mension. Based on this assessment, we decided to focus 
on the non-invariant form in this work. 

Should one choose to implement the Galilean invariant 
form instead, as in Ref. [4], there are two technical de- 



tails that should not be overlooked. First, in Eq. ( 16 ) the 
function A{p^) becomes A(p^—q^) and therefore needs to 
be evaluated outside the single-particle momentum lat- 
tice where we originally defined it (c.f. Eq. (13l). The 



extension to the larger domain should be done respecting 
the periodicity of the momentum lattice. 

Second, if such a periodicity is to be reconciled with 
the continuum form of the non-interacting dispersion re- 
lation, which we assumed io he E = p^ I2ra, we should 
impose a spherically symmetric cutoff A in momentum 
space. This last condition should have little impact on 
the results in practice, as long as the systems under con- 
sideration are somewhat dilute. From the point of view of 
computational performance, however, such a restriction 
on phase space results in important gains, particularly 
at finite temperature, where all the elements in the basis 
are evolved in imaginary time. 

Regardless of the form of the action, it is possible to 
evaluate the degree to which Galilean invariance is bro- 
ken by taking the improved transfer matrix and com- 
puting its eigenvalues in a moving frame. As long as 
Galilean invariance is respected, inserting the eigenval- 
ues into the Rummukainen-Gottlieb formula [2 5) should 
yield the same scattering phase shift as Liischer's for- 
mula. This will of course include breaking effects coming 
from both the action and the lattice itself. 

Here, we limit ourselves to comparing the eigenvalues 
obtained from the invariant and non-invariant formula- 
tions, in various frames. In Fig. [7] we show a plot of 
the relative eigenvalue difference |Ai?„/i?„| between the 
Galilean invariant and non-invariant improved transfer 
matrices, evaluated in various moving frames (i.e. non- 
zero center-of-mass momentum). As can be appreciated 
in that plot, the difference between the invariant and non- 
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invariant improvements is extremely small, which shows 
that one can use the non-invariant version of the improve- 
ment program without introducing a noticeable system- 
atic effect. 



FIG. 7. (Color online) Logarithmic plot of the relative eigen- 
value difference | A_B„/£'„| between the Galilean invariant and 
non-invariant improved transfer matrices, as a function of the 
eigenvalue index n, in various moving frames. The quantity 
Pqm denotes the center-of-mass momentum, in units oi2n/L. 
This data set corresponds to N^ = 18, t = 0.05, Nq — 3. 
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